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CONVERSION  FACTORS,  NON-SI  TO  SI  (METRIC) 
UNITS  OF  MEASUREMENT 


Non-SI  units  of  measurement  used  in  this  report  can  be  converted  to  SI  (met¬ 
ric)  units  as  follows: 


Multiply _ 

Fahrenheit  degrees 

gallons 

inches 


_ 2x_ 

5/9 

3.785412 

2.54 


_ To  Obtain _ 

Celsius  degrees  or  kelvins* 
cubic  decimetres 
centimetres 


*  To  obtain  Celsius  (C)  temperature  readings  from  Fahrenheit  (F)  readings, 
use  the  following  formula:  C  ■  (5/9) (F— 32) .  To  obtain  Kelvin  (K) 
readings,  use  K  =■  (5/9)  (F-32)  +  273.15. 
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GEOTECHNICAL  APPLICATIONS  OF  THE  SELF  POTENTIAL  (SP)  METHOD 

Report  4 

NUMERICAL  MODELING  OF  SP  ANOMALIES:  DOCUMENTATION 
OF  PROGRAM  SPPC  AND  APPLICATIONS 


PART  I :  INTRODUCTION 


Background 


1.  This  report  is  the  fourth  in  a  series  of  reports  on  geotechnical 
applications  of  the  self  potential  (SP)  method.  The  following  is  a  brief 
description  of  the  previous  reports  in  this  series: 

Report  1:  Introduces  and  describes  a  novel  technique  using  the  SP 
method  for  detecting  and  mapping  subsurface  flow  in  and 
around  sinkholes  in  karst  areas  (Erchul  1988); 

Report  2:  Extends  the  concepts  introduced  in  Report  1  to  more  general 
ground-water  flow  mapping  in  karst  areas;  investigates  the 
effects  of  environmental  variables,  such  as  rainfall, 
temperature,  soil  type,  etc.,  on  SP  surveys;  examines  the 
requirements  for  and  the  utility  of  an  automated  SP  data  col¬ 
lection  system  (Erchul  and  Slifer  1989); 

Report  3:  Investigates  the  SP  data  acquisition  requirements  for 
subsequent  quantitative  modeling  and  interpretation, 
electrode  types,  responses  and  comparisons  and  effects  of 
environmental  variables;  presents  a  recommended  field  data 
acquisition,  procedure;  presents  analytical  interpretation 
procedures  using  geometric  models;  introduces  the  interpre¬ 
tation  of  SP  data  using  numerical  modeling  and  applies  the 
method  to  modeling  the  SP  anomaly  caused  by  flow  through  an 
earth  dam  model;  presents  an  extensive,  annotated  bibliogra¬ 
phy  to  the  SP  literature  (Corwin  1989)  ; 

This  report,  Report  4,  documents  a  microcomputer  program  for  the  numerical 
modeling  of  SP  anomalies  caused  by  flow  of  fluid  or  heat  (primary  flows) 
through  a  porous  medium.  The  program  is  applied  to  the  modeling  of  SP  data 
acquired  at  Beaver  Dam,  Arkansas.  This  model  calculation  is  the  first  known 
numerical  modeling  simulation  of  SP  phenomena  caused  by  anomalous  seepage 
through  an  earth  dam  and  foundation. 
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Scope 


2.  Program  SPPC  is  a  FORTRAN  code  used  to  calculate  self-potential  (SP) 
anomalies  due  to  thermal  and  pressure  sources.  This  code  is  a  version  of  pro¬ 
gram  SPXCPL  (Sill  1983;  Sill  and  Killpack  1982)  that  has  been  modified  for  use 
on  an  IBM  compatible  personal  computer  (PC)  running  MSDOS.  It  is  a  2-1/2  d 
solution — that  is,  the  program  restricts  the  user  to  a  two-dimensional  distri¬ 
bution  of  physical  properties  but  allows  for  a  three-dimensional  source 
function. 

3.  The  program  calculates  voltages  and  primary  potentials  (pressure  or 
temperature)  throughout  the  two-dimensional  space.  The  final  solution  is 
generated  using  a  three  part  process.  First,  the  pressure  or  temperature 
distribution  is  determined  from  the  source  parameters  and  the  permeability  or 
thermal  conductivity  distribution.  After  this  primary  potential  (pressure  or 
temperature)  has  been  calculated,  the  locations  and  strengths  of  the  induced 
electrical  sources  are  determined.  These  electrical  sources  are  determined 
from  the  primary  potentials  (pressure  or  temperature)  and  the  distribution  of 
the  cross-coupling  coefficients  which  relate  the  electrical  and  the  hydraulic 
or  thermal  potentials.  At  this  point  the  problem  is  reduced  to  electrical 
conduction  equation.  That  is,  the  final  SP  voltages  are  computed  from  the 
induced  electrical  sources  and  the  two-dimensional  resistivity  distribution. 

4.  Although  the  code  is  fairly  simple  to  operate,  it  is  not  quite  so 
simple  to  use  it  to  model  physical  problems.  The  difficulty  comes  from  the 
necessity  to  know  two-dimensional  distributions  of  primary  impedance  (perme¬ 
ability  or  thermal  conductivity),  cross-coupling  coefficients,  and  electrical 
resistivity.  Some  practice  is  also  required  to  assign  appropriate  boundary 
conditions  and  source  parameters  and  to  determine  if  the  computed  solution  is 
physically  reasonable.  This  manual  will  guide  the  user  through  several 
examples  so  that  the  modeling  process  is  clear. 

3.  The  manual  is  divided  into  three  sections  in  addition  to  Appendix  A. 
The  first  section  briefly  develops  the  field  equations  and  describes  the 
method  of  solution.  It  should  provide  the  user  with  a  basic  understanding  of 
the  workings  of  the  program.  More  complete  treatments  of  these  topics  are 
given  in  Sill  (1983)  and  Marshall  and  Madden  (1959).  The  second  section 
describes  the  input-output  operations  for  running  the  program  and  provides 
some  details  on  the  actual  code.  A  sample  problem  is  used  to  show  how  a 
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physical  model  may  be  translated  into  the  appropriate  input  file  and  how  the 
output  can  be  examined  to  see  if  it  is  producing  a  physically  reasonable 
result.  The  final  section  of  the  manual  presents  a  field  example  case  for  the 
modeling  of  SP  associated  with  a  leaky  dam.  This  section  should  provide  the 
user  with  a  glimpse  at  the  power  and  limitations  of  this  code.  Appendix  A 
contains  some  values  for  permeability,  thermal  conductivity,  thermal  and 
hydraulic  cross-coupling  coefficients  and  electrical  resistivity  of  earth 
materials  that  might  be  reasonably  encountered  in  field  surveys.  The  tables 
are  brief  because,  for  many  rocks,  some  of  these  physical  property  values  are 
not  well  known. 

6.  Throughout  this  manual  it  is  assumed  that  the  user  has  a  working 
knowledge  of  the  operation  of  an  IBM  PC  compatible  computer.  The  program 
requires  the  user  to  modify  and  store  input  files  and  to  retrieve  and  print  or 
plot  output  files.  Included  with  this  manual  is  a  5-1/4-in.*  disc  containing 
five  files.  These  include  source  and  executable  code,  and  sample  input, 
output  and  plotting  files.  We  strongly  recommend  that  you  make  a  copy  of  this 
disc  when  you  first  begin  to  use  the  program  and  use  the  original  disc  as  a 
back-up.  There  are  no  restrictions  on  modifying  the  program. 


*  A  table  of  factors  for  converting  non-SI  to  SI  (metric)  units  of  measure¬ 
ment  is  presented  on  page  3. 


6 


PART  II:  THEORETICAL  BACKGROUND 


7.  The  basic  principle  relating  self  potential  and  fluid  or  heat  flow 
processes  is  that  the  flow  of  electrical  current  and  heat  or  fluid  are 
coupled.  That  is,  there  are  electrical  potentials  <|>  related  to  fluid  or 
heat  flow  processes  and  hydraulic  or  thermal  P  associated  with  electrical 
current  flow.  The  mathematical  relationship  between  the  potentials  is  given  in 
Equations  (1),  (2)  (Nourbehecht  1963;  Marshall  and  Madden  1959). 


Q  =  -  CUVP  -  C12V* 


(1) 


J  -  -  Cn7P  -  C22V<j>  (2) 

where 

Q  and  J  =  primary  (heat  or  fluid)  flow  and  current  flow  vectors, 
respectively 

Cj  =  K,  the  hydraulic  or  thermal  conductivity 

^22  =  C’  t^e  e^-ectrica^  conductivity 

Cf2  =  C^,  the  cross-coupling  coefficients 

8.  The  cross-coupling  terms  in  Equations  (1),  (2)  are  typically  much 
smaller  than  the  primary  flow  terms.  Fluid  or  heat  flow  processes  do  not 
generate  very  large  electrical  currents  and  imposition  of  electrical  potential 
does  not  generate  large  fluid  flows  (although  this  method  which  is  called 
electroosmosis  has  been  used  for  dewatering  low  permeability  material).  Note 
that  we  neglect  the  cross-coupling  terms  in  Equations  (1),  (2)  the  equations 
decouple  into  the  more  familiar  Darcy's  law  and  Ohm's  law. 

9.  In  the  absence  of  external  electrical  current  sources  we  note  that 
the  current  is  divergenceless 


V  •  J  '  ^  «  0 

dt 


where  q  is  electrical  charge.  Applying  this  relation  to  Equation  2  in 
addition  to  some  algebra,  we  obtain 
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(4) 


-  V  X  (0V<j>)  =  VC  x  Vp+  C21V2P 

10.  This  is  the  classical  dc  conduction  equation  with  electrical 

potential  and  conductivity  on  the  left-hand  side  and  the  induced  current 

sources  on  the  right-hand  side.  Equation  (4)  shows  that  induced  current 

sources  occur  if  (a)  the  cross-coupling  coefficient  C_.  changes  in  the 

^  2 

direction  of  primary  flow  and/or  (b)  there  are  primary  flow  sources  (V  p  ^  0) . 

11.  Because  the  effects  of  secondary  electrical  potential  on  fluid  or 
heat  flow  are  small,  the  pressure  (or  temperature)  distribution  may  be 
obtained  by  a  straightforward  application  of  Darcy's  law  or  the  heat  conduc¬ 
tion  equation.  This  primary  potential  (pressure  or  temperature)  is  then  used 
in  combination  with  the  coupling  coefficient  distribution  to  calculate  the 
electrical  source  terms  on  the  right-hand  side  of  Equation  (4).  Once  the 
source  terms  are  known,  the  electrical  potential  may  be  calculated  from  Equa¬ 
tion  (4)using  numerical  methods. 

Numerical  Modeling  Scheme 

12.  Program  SPPC  uses  a  transmission  surface  analogy  to  solve  Equa¬ 
tion  (4).  The  transmission  surface  approach  draws  an  analogy  between  the  physi¬ 
cal  process  to  be  modeled  and  a  power  transmission  network  (Swift  1967;  Madden 
1971).  With  this  approach  the  model  space  is  transformed  into  a  lumped  ele¬ 
ment  network  and  the  terms  of  Equation  (4),  for  example,  are  replaced  by 
impedances  and  admittances  at  each  node.  By  imposing  Kirchoff's  circuit  law 
and  the  continuity  of  current  at  each  node  of  the  network,  we  end  up  with  an 
array  of  finite  difference  equations. 

13.  Because  the  program  uses  a  three-dimensional  source  function  and 
two-dimensional  physical  property  distributions,  the  equations  must  first  be 
Fourier  transformed  to  eliminate  the  three-dimensional  source  dependence 
(Swift  1967;  Sill  1983).  The  two-dimensional  finite  difference  equations  (in 
Fourier  domain)  are  then  solved  implicitly  using  Greenfield's  algorithm  (Swift 
1967).  Solutions  are  obtained  at  several  Fourier  wavenumbers  and  the  solution 
in  three-dimensional  space  is  obtained  by  inverse  Fourier  transform. 
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PART  III:  PROGRAM  OPERATIONS 


File  Preparation  and  a  Sample  Problem 


Problem  description 

14.  In  this  section,  program  SPPC  is  used  to  calculate  the  self  poten¬ 
tial  anomaly  for  a  simple  problem  where  an  analytical  solution  is  available. 
This  provides  the  user  with  an  example  of  how  the  model  parameters  are  trans¬ 
lated  into  a  computer  input  deck.  Comparing  the  analytical  and  numerical 
solutions  for  this  problem  also  provides  a  validation  of  the  computer  code. 

15.  The  problem  considered  is  a  point  pressure  source  in  a  homogeneous 
half space.  For  the  case  the  input  deck  is  completely  described  and  the  model 
parameters  are  matched  one  by  one  to  the  program  variables  in  the  input  deck. 
The  output  is  then  compared  with  the  analytical  results. 

16.  Figure  1  schematically  illustrates  the  problem  of  a  point  pressure 
source  in  a  homogeneous  half space.  This  may  correspond  to  calculating  the  SP 
anomaly  around  an  injection  well,  for  example.  An  analytical  solution  of  this 


image  wall 

O 

r 

z 


observation  point 
_ 


z 

A 

true  well 


Figure  1 .  Schematic  model  of  a  pressure  source  in  a 
homogeneous  halfspace 
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problem  is  easily  obtained  by  placing  an  image  source  above  the  surface  at  a 
height  equal  to  the  depth  of  the  true  source.  The  superposition  of  the  true 
source  and  the  image  source  satisfies  the  governing  equation  (Darcy's  law)  and 
also  matches  the  boundary  condition  at  the  surface  (no  vertical  flow).  The 
derivation  below  follows  the  treatment  of  Sill  and  Johng  (1981). 

Analytical  solution 

17.  The  pressure  distribution  may  be  obtained  by  applying  Darcy's  law 
to  the  true  source  and  the  image  source 

Q  “  ~  VP  (5) 


where 

Q  =  flow  rate 

k  =  intrinsic  permeability,* 

A  =  cross-sectional  area 
p  *  fluid  viscosity 
P  =  pressure 

For  a  single  point  source  in  a  homogeneous  earth,  Equation  (5)  may  be  solved 
for  the  pressure  to  give 


P(r) 


=  _2h_ 

4Trrk 


(6) 


where  r  is  the  distance  from  the  source  to  the  observation  point. 

18.  For  a  homogeneous  half space,  the  pressure  from  the  image  well  must 
also  be  calculated. 


P(r.r') 


Jk_  +  Qp 

4irrk  4iTr'k 


(7) 


*  The  hydraulic  conductivity  K  is  related  to  the  intrinsic  permeability  k 
by  K  =  kog/p  ,  where  g  is  the  acceleration  of  gravity  and  a  and  p 
are  the  density  and  viscosity,  respectively  of  the  fluid. 
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At  the  surface  r  and  r'  are  equal  In  magnitude  and,  therefore,  the  solu¬ 
tion  at  the  surface  is 


where 


(8) 


x  «=  horizontal  distance  from  the  source  to  the  observation  point 
Zq  =*  source  depth 


19.  The  electrical  potential  at  the  surface  follows  from  Equation  (2). 
This  also  has  a  simple  form 


where 


(9) 


p  =  resistivity  of  the  halfspace 
Cq  -  cross-coupling  coefficient 

20.  One  difficulty  in  SP  calculations  is  maintaining  consistency  in  the 
units.  Because  fluid  flow  units  are  usually  expressed  in  the  cgs  system 
(i.e.,  darcies  and  centipoise  for  k  and  p)  and  electrical  units  are  most  often 
given  in  the  MKS  system  (i.e,  ohm-metres  and  volts  for  p  and  V)  unit  conver¬ 
sions  are  needed  to  obtain  results  in  one  system  or  the  other.  Program  SPPC 
utilizes  the  standard  set  of  units  for  fluid  or  heat  flow  and  electrical 
parameters  and  makes  unit  conversions  within  the  program.  The  units  for  the 
program  input  are  given  in  Table  1.  Remember  that  all  input  values  must  be 
given  in  these  units  or  the  solution  is  invalid. 

Program  input  and  output 

21.  Figure  2  illustrates  the  input  deck  for  program  SPPC  for  the 
problem  of  a  point  fluid  source  of  1  i/s  in  a  halfspace  at  a  depth  of  22.5  m. 
The  halfspace  parameters  are  given  as  follows:  p  ■  10  ohm-m,  -  mv/atm, 

k  ■  10.0  millidarcies .  These  are  typical  values  for  unconsolidated  sedimen¬ 
tary  rocks  (Appendix  A).  The  section  below  describes  the  input  parameters  in 
detail. 
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Table  1 

Units  In  Program  SPPC 


Units  for  Sources 


Fluid  flow 
Heat  flow 


liters/sec 

watt 

Units  of  Physical  Properties 


Permeability* 

Cross-coupling  coefficient  (fluid) 
Cross-coupling  coefficient  (thermal) 
Resistivity 


_q  2 

darcies  (1  darcy  **  9.87  x  10  cm  ) 
mvolt/ atmosphere 
mv/°  C 
ohm-m 


Calculated  Data 


Pressure 
Temperature 
Electrical  potential 


atmosphere  (1  atm 

0  C 

volt 


=  1.013  -  x  10  Pa) 


*  For  water  at  20°  C,  1  darcy  (permeability  unit)  is  equivalent  to 
_4 


9.613  x  io  cm/s  (hydraulic  conductivity  unit). 


HALFSPACE 

112 

10.0 

1 

26 

6 

1.0 

0 

PHYSICAL  PROPERTIES 


!  PROFILE  TITLE 

!  1=  PRESS  2 “THERM,  LINE  SOURCE  (0  OR  1),  LLENGTH 

!  LENGTH  9m)  OF  MODEL  UNIT 

!  #  OF  SOURCES 

!  X  LOCATION  OF  SOURCES 

!  Z  LOCATION  OF  SOURCES 

!  STRENGTH  OF  SOURCES  (Li/S  OR  WATTS) 

!  CHANGE  ARRAY  WEIGHTING  Z  (0  OR  1) 


5 

PERM/TC 

CCP 

RESIS 

1 

10.0E-3 

20.0 

10.0 

1 

PARAMETER  VALUES 

FOR  UNITS 

2 

10.0E-3 

20.0 

10.0 

J 

PERM 

(DARCIES) 

3 

10.0E-3 

20.0 

10.0 

1 

CCP 

(MV/ ATM) 

4 

10.0E-3 

20.0 

10.0 

1 

RESIS 

(OHM-M) 

5 

10.0E-3 

20.0 

10.0 

5  10 

15  20 

25  30 

35 

40 

45  50 

-- 

- + - + — 

— + - +— 

- +C - +— 

- +. 

- + — 

!  Z-SPACING 

11111111111111111111111111111111111111111111111111111 

21111111111111111111111111111111111111111111111111111 

31111111111111111111111111111111111111111111111111111 

41111111111111111111111111111111111111111111111111111 

51111111111111111111111111111111111111111111111111111 

6111111111111111111111111]G)llllllllllllllllllllllllll 

71111111111111111111111111111111111111111111111111111 

81111111111111111111111111111111111111111111111111111 

91111111111111111111111111111111111111111111111111111 

101111111111111111111111111111111111111111111111111111 

111111111111111111111111111111111111111111111111111111 


.25 

.25 

.25 

.50 

.50 

.50 

.50 

1.0 

1.0 

1.0 

2.0 


Figure  2.  Sample  SPPC  input  file  for  the  case  of  a  point  pressure 
source  in  a  homogeneous  halfspace 
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Input  deck 

22.  Program  SPPC  uses  a  mesh-type  input  which  is  common  to  many  numeri¬ 
cal  modeling  programs.  It  offers  great  flexibility  but  requires  some  adjust¬ 
ment  if  you  are  unfamiliar  with  this  type  of  input.  With  this  program  the 
locations  of  all  source  and  measurement  points  are  referred  to  coordinates  in 
the  input  mesh.  You  must  therefore  keep  track  of  the  "true"  locations  by 
using  a  distance  scale  factor.  For  example,  if  the  scale  factor  is  100  m, 
then  each  "unit"  within  the  mesh  represents  100  m  on  or  within  the  earth.  A 
scale  "unit"  corresponds  to  4  mesh  points  in  the  X  direction  and  between  0.5 
and  0.4  points  in  the  Z  direction  depending  on  the  depth. 

23.  The  program  works  by  assigning  each  individual  mesh  coordinate  a 
set  of  physical  properties:  permeability  of  thermal  conductivity,  cross¬ 
coupling  coefficient,  and  electrical  resistivity.  It  assumes  that  these 
parameters  are  constant  within  the  mesh  block.  The  model  therefore  consists 
of  a  two-dimensional  array  of  these  mesh-blocks.  If  the  blocks  are  small 
enough,  then  the  discontinuous  model  will  approximate  the  more  continuous 
earth,  but  keep  in  mind  that  the  voltages,  temperature,  and  pressures  are 
known  only  at  these  mesh  points. 

24.  The  Input  is  format  free,  that  is  input  values  do  not  have  to 
appear  in  specific  columns  in  the  deck.  Comments  appearing  in  the  input  deck 
are  preceded  by!.  Note  that  the  two  lines  giving  the  linear  scale  must  be 
included  in  the  input  deck  as  the  program  expects  to  see  them.  The  individual 
input  variables  are  described  below: 
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LINE  1  :  TITLE 


Profile  title  (15  characters  max) 


This  input  consists  of  alphanumeric  characters  corresponding  to  the  title  of  the  profile.  The  program 
will  ignore  characters  beyond  the  maximum  limit  of  15. 

LINE  2  :ITYPEJLINEXLENGTH  Type  Sag,  line  source  flag.line  length 

ITYPE=1  for  pressures  sources, =2  for  thermal  sources.  Values  other  than  1  or  2  are  reset  by  the  pro¬ 
gram  to  1  (pressure  sources). 

QJNE=0  Calculate  for  point  sources  only,  =1  Calculate  line  sources  also. 

I  .LENGTH  *  Length  of  the  line  source  in  model  units;  the  source  is  centered  around  Y=0.  The  source 
is  uniformly  distributed  throughout  the  line.  NOTE;  The  output  for  a  line  source  is  given  in  cross- 
section  form  at  the  end  of  the  file.  The  cross-sections  and  areal  plots  at  the  front  of  the  output  refer  to 
point  sources  only. 

LINE  3  rASCALE  Length  in  meters  of  one  model  unit 

Note  that  1  model  unit  corresponds  to  4  mesh  points  in  the  X  direction.  In  the  Z  direction  the  mesh 
spacing  depends  on  the  Z  coordinate.  For  the  default  spacing,  given  in  Figure  2,  the  first  3  Z  coordi¬ 
nates  are  spaced  .25  units  the  next  4  are  spaced  .5  units,  the  next  3  are  spaced  1  unit  and  the  final  spac¬ 
ing  is  2  units. 

LINE  4  :NT  Number  of  sources  (pressure  or  thermal) 

The  maximum  number  of  sources  allowed  is  9. 

LINE  5  :NTX(I),  1=1, NT  X  coordinate  of  each  source 

The  X  location  of  each  source  is  given  as  a  coordinate  in  the  mesh  (see  sample  deck  in  Figure  2). 

LINE  6  :NTZ(I),  1=1, NT  Z  coordinate  of  each  source 

The  Z  location  of  each  source  is  given  as  a  coordinate  in  the  mesh  (see  sample  deck  in  Figure  2). 

LINE  7  ;  TS(D,  1=1, NT  Strengths  of  the  various  sources 
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Source  strengths  for  pressure  sources  are  given  in  liters  per  sec.  (Note  that  one  liter  per  sec  is  roughly 
equal  to  16  gallons  per  minute).  The  source  strength  for  thermal  sources  is  given  in  watts. 


LINE  8:  IZ  Flag  for  changing  Z  mesh 

If  IZ=1  then  a  new  mesh  is  input  for  Z.  if  IZ=0  then  the  default  mesh  is  used.  The  default  mesh  is 
shown  in  Figures  2  and  3. 

LINE  9:DELV(T)J=1,11  New  mesh  spacing  for  Z 

If  12=1  then  input  11  new  Z  mesh  spacings.  Spacings  do  not  need  to  increase  with  depth  but  the  final 
spacings  need  to  be  large  to  insure  an  accurate  solution. 

LINE  10:NU  Number  of  individual  model  units 

The  maximum  allowable  number  of  units  is  9.  In  the  example  shown  in  Figure  2  NU=5. 

LINE  11,I1+NU:IJPERM(I),CCP(I),RES(I)  Physical  property  values 

PERM(I)  permeability.  CCP(I)  cross-coupling  coefficient  and,  RES(I)  resistivity  of  the  various  model 

units. 

PERM(I),  is  the  permeability,  in  darcies.  or  thermal  conductivity  in  mwatts/  m-deg  C. 

CCP(D.  is  the  cross-coupling  coefficient  of  the  various  units  in  mv/atm  or  mvolt/  deg  C. 

R£S(I),  is  the  resistivity  of  the  various  units  in  ohm-meters. 

In  the  example  shown  in  Figure  2  all  five  units  have  the  same  permeability,  cross-coupling  coefficient 
and  resistivity. 

LINES  I2+NU2J+NU  Gridded  mesh 

Each  number  corresponds  to  mesh  unit  which  has  the  three  physical  property  values  as  shown  above. 
This  mesh  is  the  repesentation  of  the  input  model,  the  center  of  the  mesh  (node  26)  is  given  as  the 
coordinate  0  in  the  output.  In  the  X-direcuon  all  mesh  units  are  spaced  .25  units  apart  in  the  Z  direction 
the  mesh  spacing  vanes  from  .25  units  near  the  surface  to  2  units  at  depth. 
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Running  the  program 

25.  After  preparing  an  input  deck  as  shown  above  and  saving  it  on  the 
disk  of  your  PC  you  may  run  the  program  simply  by  typing  "SPPC  fn"  where  "fn" 
designates  the  name  of  the  input  file  to  be  computed  (the  .dat  extension  is 
automatically  provided).  To  run  the  sample  problem  provided  simply  type  SPPC 
HALF.  The  program  will  automatically  assign  HALF.DAT  as  the  input  file  name, 
HALF. RES  as  the  output  file  name  and  HALF.PLT  as  the  plotting  file  name.  For 
this  sample  problem  the  program  requires  about  8  min  to  run  on  an  IBM-PCAT; 
additional  time  is  required  if  more  than  one  source  is  used. 

26.  The  plotting  file  contains  only  the  model  name,  the  distance  scale 
factor,  and  the  primary  potential  and  voltage  at  individual  mesh  points  on  the 
surface.  No  plotting  software  is  provided  with  this  manual  but  a  number  of 
commercial  plotting  programs  written  for  PC  compatible  microcomputers  accept 
this  type  of  file  for  input. 

27.  The  output  for  the  above  example  is  shown  in  Figure  3.  The  first 
section  of  the  output  file  contains  a  summary  of  the  input  parameters  and  a 
cross  section  of  the  physical  parameter  distribution.  The  results  of  the 
model  calculation  are  then  given  as  a  profile  of  primary  potentials  and  volt¬ 
ages  at  the  surface  and  then  as  detailed  profiles  of  primary  potentials  and 
voltages.  The  third  section  is  an  areal  distribution  (plan  view)  of  the  sur¬ 
face  voltage.  The  final  section  of  the  output  gives  cross-sections  of  the 
voltages  from  a  line  source.  For  these  sections  the  line  sources  are  centered 
about  Y«=0  and  the  strength  of  the  source  in  distributed  evenly  along  the 
length  of  the  line.  The  cross  sections  will  appear  only  if  ILINE  =  1  in 

line  1  of  the  input  file. 

28.  A  comparison  of  the  numerical  model  and  the  analytical  results  for 
points  at  the  surface  is  shown  in  Figure  4.  Examining  both  the  pressure  and 
voltage  curves  in  Figure  4  we  find  a  good  agreement  between  numerical  and 
theoretical  results.  Errors  of  1  to  5  percent  are  observed  at  the  farther 
mesh  point.  These  are  acceptable  for  this  type  of  numerical  modeling  although 
they  may  be  reduced  by  providing  a  finer  mesh  discretization. 

29.  Before  accepting  the  numerical  solution  it  is  important  to  deter¬ 
mine  that  the  output  is  providing  physically  reasonable  results.  Note,  for 
example,  in  the  cross-sectional  plots,  the  pressure  is  highest  near  the  source 
location.  Also  notice,  from  the  surface  areal  output,  that  the  potential  is 
symmetrical  about  the  surface  coordinate  of  the  source  point. 
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Figure  3.  SPPC  output  file  for  a  point  pressure  source  in  a  homogeneous  halfspace  (Continued) 


Figure  3.  (Concluded) 


PRESSURE  (atm) 


Program  Notes 


30.  A  few  rules  of  thumb  are  useful  in  applying  this  program  to  model 
physical  problems.  First,  when  constructing  an  input  mesh,  it  is  a  good  prac¬ 
tice  to  sketch  out  your  model  with  paper  and  pencil  before  you  type  in  the 
input.  This  allows  a  check  on  the  input  to  see  if  it  is  consistent  with  the 
model  that  you  want.  Second,  when  constructing  an  input  mesh,  it  is  recom¬ 
mended  that  you  begin  at  the  center  of  the  model  and  work  outwards.  The 
center  point  is  located  at  node  26  and  remember  that  the  X  nodes  are  spaced 
0.25  units  apart  throughout  the  visible  input  mesh.  Although  the  complete  X 
mesh  actually  consists  of  68  nodes,  only  the  central  52  nodes  are  used  to 
input  models;  the  remaining  nodes  extend  the  mesh  laterally  to  eliminate 
numerical  errors.  The  output  file  provides  an  expanded  plot  of  the  voltage 
and  primary  potential  for  the  central  part  of  the  mesh;  therefore,  it  is  use¬ 
ful  to  have  your  model  centered  around  this  region.  The  expanded  region  is 
located  between  -2  and  2  units  and  the  voltages  and  primary  potentials  are 
plotted  at  each  node  (0.25  unit)  in  this  region;  whereas,  outside  of  this 
area,  they  are  plotted  every  4  nodes  (1  unit). 

31.  Modeling  SP  data  for  a  thermal  source  is  often  more  difficult  than 
fluid  flow  modeling  because  the  parameters  are  poorly  constrained.  Because 
the  nature  of  heat  sources  are  usually  less  well  known  than  fluid  flow 
sources,  it  is  sometimes  useful  to  initially  assume  a  source  strength 
(100,000  watts  for  example)  and  examine  the  temperatures  in  the  output  file  to 
see  if  they  are  reasonable  for  your  model.  See  the  papers  by  Sill  (1983)  and 
Sill  and  Johng  (1981)  for  examples  of  thermal  source  modeling. 
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PART  IV:  FIELD  EXAMPLE 


32.  The  section  below  illustrates  the  use  of  this  code  for  a  more  com¬ 
plex  example.  In  this  case  the  code  is  used  to  develop  a  model  for  a  leaky 
damsite  where  a  field  SP  survey  has  been  conducted.  This  example  will  high¬ 
light  some  of  the  considerations  in  using  the  code  to  model  field  results. 

33.  Figure  5  shows  a  self  potential  anomaly  map  for  the  Beaver  damsite 
in  the  state  of  Arkansas  in  the  southern  part  of  the  USA  (Llopis  and  Butler 
1988).  This  is  an  earthen  dam  through  which  significant  seepage  occurs  in  the 
region  labeled  "new  wet  area".  The  seepage  in  the  southern  part  of  the  survey 
area  is  thought  to  be  controlled  by  an  east-west  trending  fault  zone  that 
obliquely  crosses  the  dam  structure.  The  fault  forms  the  southern  boundary  of 
a  graben  structure  and  is  associated  with  a  string  of  negative  SP  anomalies 
(Butler  1988). 

34.  The  SP  map  shows  some  of  tha  general  characteristics  of  a  leakage 
problem  (Figure  5).  That  is,  there  are  negative  anomalies  where  fluid  is 
leaking  into  the  dam  and  reservoir  boundary  (near  A)  and  positive  anomalies 
over  the  surface  discharge  areas  (new  wet  area).  The  association  of  the 
anomalies  with  the  known  fault  zone  strongly  suggests  that  the  fault  provides 
a  low  impedance  leakage  pathway  for  fluids.  The  map  also  shows  several  other 
anomalies,  some  of  which  are  related  to  other  leakages  and  some  are  due  to 
topography  which  can  have  a  significant  effect  (Corwin  1988). 

35.  We  selected  a  single  profile  from  the  map  (A-A')  that  connects 
source  and  seep  areas  and  fit  the  observed  data  to  calculate  data  for  a  two- 
dimensional  geometry.  The  permeability  and  resistivity  distributions  for  the 
model  shown  in  Figure  6  were  obtained  from  field  measurements  (Butler  1988). 
Values  given  for  the  cross-coupling  coefficients  are  "educated  guesses"  based 
on  measurements  reported  in  Ishido  and  Mizutani  (1981)  and  Noubehecht  (1963) 
(Appendix  A) .  Figure  6  is  a  very  simple  model  for  the  dam  and  only  crudely 
matches  the  geometry  and  construction  details  of  the  dike.  The  model  features 
single  source  and  seep  locations  and  a  single  coupling-coefficient  contrast 
corresponding  to  where  the  graben  fault  crosses  the  profile  line.  The  magni¬ 
tude  of  the  leak  is  about  0.5  £/ sec  (about  8  gal/min)  and  about  one-half  of 
this  amount  is  discharging  in  the  new  wet  area. 

36.  The  input  deck  corresponding  to  the  Beaver  Dam  model  is  given  in 
Figure  7;  it  has  several  interesting  features.  Since  water  must  flow  out  of 
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Figure  7.  Input  file  for  the  Beaver  Dam  field  example 


the  dam  in  the  downstream  direction,  "permeability"  of  the  water  in  the 
reservoir  (property  set  2)  is  assumed  to  be  very  low.  This  ensures  that  water 
from  the  source  will  not  flow  back  into  the  reservoir.  Also,  since  the  dam  is 
a  topographic  structure,  measurements  must  be  made  within  the  mesh  and  not  on 
the  surface  (see  marked  input  in  Figure  7). 

37.  For  the  Beaver  Dam  model  15  iterations  of  forward  modeling  were 
required  to  match  field  results.  The  fit  between  calculated  and  observed  SP 
is  shown  in  Figure  8.  A  very  good  match  was  achieved  even  though  the  assumed 
model  is  considerably  simpler  than  the  known  geology.  This  suggests  that  the 
SP  anomalies  are  dominantly  controlled  by  the  locations  and  magnitudes  of  the 
sources  and  seeps  and  not  by  complex  flow  geometry  or  rock-type  variations. 
However,  this  is  not  a  unique  model  and  the  parameters  are  not  well 
constrained. 
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PART  V:  CONCLUSIONS  AND  RECOMMENDATIONS 


38.  This  paper  describes  the  theory  and  operation  of  computer  program 
SPPC.  The  program  calculates  SP  anomalies  due  to  fluid  and  heat  flow 
processes  in  a  two-dimensional  earth.  The  model  uses  a  mesh-type  input  and 
may  require  some  practice  before  modeling  can  be  done  routinely;  it  runs  on  an 
IBM  PC-AT  computer  in  about  8  min. 

39.  This  program  is  potentially  a  powerful  tool  for  geoscientists  and 
engineers  using  the  SP  method  to  locate  areas  of  subsurface  ground-water  flow 
such  as  leaks  in  earthen  dams,  dikes,  or  evaporation  ponds.  It  is  also  a 
useful  tool  in  geothermal  exploration.  Future  applications  may  include  such 
diverse  problems  as  nuclear  repository  monitoring  or  tracking  of  subsurface 
contaminant  flow.  Note  that  the  program,  in  its  present  form,  is  designed  for 
use  on  simple  problems.  Due  to  the  memory  restrictions  and  limited  computing 
power  on  a  PC,  the  mesh  is  quite  small  so  only  fairly  simple  models  are 
possible.  Also,  the  potential  user  is  encouraged  to  consult  the  References 
and  particulary  Report  3  (Corwin  1989)  for  guidance  on  obtaining  required 
model  input  parameters. 
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APPENDIX  A:  PHYSICAL  PARAMETERS 


1 .  Some  representative  physical  property  values  that  may  be  used  in 
modeling  SP  data  are  provided  in  Table  A1 .  Included  are  values  for  resistiv¬ 
ity,  permeability,  thermal  conductivity,  electrokinetic  (fluid  flow),  and 
thermoelectric  cross-coupling  parameters.  These  data  were  compiled  from  a 
number  of  different  sources;  the  values  given  are  given  as  broad  ranges. 

Values  for  permeability  are  compiled  from  Freeze  and  Cherry  (1979)  thermal 
conductivity  information  is  obtained  in  Nourbehecht  (1963)  and  Welch  et  al. 
(1981).  Values  for  electrical  resistivity  are  obtained  from  Keller  and 
Friscbknecht  (1966),  and  Keller  (1988).  Cross-coupling  coefficients  are 
obtained  from  Nourbehecht  (1963),  Ishido  and  Mitzutani  (1981)  and  Fitterman 
and  Corwin  (1982). 

2.  The  effects  of  porosity,  pore  water  salinity,  and  specific  mineral 
composition  are  all  important  in  determining  these  physical  properties.  These 
effects  account  for  the  wide  ranges  of  values  for  these  physical  parameters. 

In  general,  the  fluid  flow  coupling  coefficients  increase  with  porosity  and 
decreases  with  pore  water  salinity.  The  resistivity  decreases  with  increased 
porosity,  pore  water  salinity,  and  temperature.  The  thermal  coupling  coeffi¬ 
cient  decreases  with  increased  porosity  and  increased  fracture  density. 

3.  The  parameters  units  are  given  as  follows:  resistivity  values  are 
given  In  ohm-n,  thermal  conductivities  are  in  watts/m-deg  C,  thermal  coupling 
coefficients  (CCth)  are  in  mv/deg  C,  permeabilities  (k)  are  in  millidarcies 
(md) ,  and  fluid  coupling  coefficient  (CCfl)  are  in  mvolts/atm.  This  table  is 
far  from  complete  and  certainly  for  many  applications  it  will  prove  to  be 
inadequate.  It  does  provide  the  user  with  a  starting  point  for  modeling  a 
number  of  physical  systems. 
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